source('00_util_scripts/mod_bplot.R')
library(clusterProfiler)
library(enrichplot)

# GSEA ---------
stts.m2hvl.deg <-
  read_csv('mission/SLE_TRPM2_MfMo/results/SLE.status.m2hvl.deg.csv')

flare.gsebp <- stts.m2hvl.deg |>
  filter(cluster == 'Flare', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL')

flare.gsebp <- flare.gsebp |>
  simplify()

flare.gsebp.tb <- flare.gsebp@result |>
  as_tibble()

flare.gsebp.tb |>
  slice_min(p.adjust, n = 9) |>
  mutate(term = fct_reorder(Description, NES)) |>
  tidyplot(x = NES, y = Description, color = p.adjust) |>
  add_mean_bar() |>
  sort_y_axis_labels() +
  scale_fill_distiller(palette = 'Reds') +
  labs(title = 'GO BP pathway enrichment of flare SLE:\nTRPM2-hi mono vs TRPM2-lo mono',
       y = 'Pathway', x = 'Normalized Enrichment Score')

publish_source_plot('SLE.pbmc.m2hvl.gogse', width = 100)

flare.gsebp.tb |>
  write_csv('mission/SLE_TRPM2_MfMo/flare.m2hvl.gsebp.csv')

flare.gsebp.tb <-
  read_csv('mission/SLE_TRPM2_MfMo/flare.m2hvl.gsebp.csv')

chemotax.gene <- flare.gsebp.tb |>
  filter(str_detect(Description, '^chemotaxis')) |>
  pull(core_enrichment) |>
  str_split_1('/')

stts.m2hvl.deg |>
  filter(gene %in% chemotax.gene, cluster == 'Flare') |>
  select(-cluster) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/taxis.pathway.gene.csv')

inflam.gene <- flare.gsebp.tb |>
  filter(str_detect(Description, '^inflamma')) |>
  pull(core_enrichment) |>
  str_split_1('/')

# TRPM2 chipseq database ---------
m2_cistrome_tf <- read_tsv('mission/SLE_TRPM2_MfMo/data/trpm2_cistrome_chipseq.tsv')

m2_cistrome_tf |>
  slice_max(score, by = c(factor, ref, cell_line)) |>
  mutate(dataset = str_c(factor, '_', ref) |> fct_reorder(score),
         Cell_origin = ifelse(cell_line %in% c('THP-1','U-937'), 'Monocytic', 'Other')) |>
  ggplot(aes(score, dataset, fill = cell_line, color = Cell_origin)) +
  geom_col(linewidth = 1) +
  scale_fill_brewer(palette = 'Paired', na.value = 'white') +
  scale_color_manual(values = c('red','#ffff'), na.value = 'white') +
  theme_classic() +
  labs(title = 'TF with regulatory potential on TRPM2 gene',
       subtitle = 'Cistrome ChIP-seq atlas')

m2_cistrome_tf |>
  slice_max(score, by = c(factor, cell_line)) |>
  mutate(factor = fct_reorder(factor, score, max)) |>
  ggplot(aes(score, factor, color = cell_line, fill = cell_line)) +
  geom_point(size = 2) +
  geom_col(width = .1, position = 'identity') +
  scale_color_brewer(palette = 'Paired') +
  scale_fill_brewer(palette = 'Paired') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'TF with regulatory potential on TRPM2 gene',
       subtitle = 'Cistrome ChIP-seq atlas', x = 'ChIP-seq binding score',
       y = 'Transcription factor') +
  theme_jpub

publish_source_plot('trpm2.chipseq.tf')

## beta output -----------
cistrome_hs_chipseq_beta <-
list.files('~/append-ssd/alaria2/beta_out/', full.names = T) |>
  map(read_tsv, id = 'file', comment ='# ', col_types = 'ciicdcc', .progress = T)

cistrome_hs_chipseq_res <-
cistrome_hs_chipseq_beta |> keep(\(x)nrow(x) > 0)

cistrome_hs_chipseq_res <-
cistrome_hs_chipseq_res |> list_rbind()

cistrome_hs_chipseq_res <- cistrome_hs_chipseq_res |>
  mutate(Externel_id = str_extract(file, 'GSM\\d+'), .keep = 'unused')

cistrome_hs_chipseq_100 <- cistrome_hs_chipseq_res |>
  slice_max(Score, n = 100, by = GeneSymbol)

cistr_hs_chipseq_meta <- read_delim('00_util_scripts/ref/hg38_tranfac_QC.txt')

cistrome_hs_chipseq_100 <- cistrome_hs_chipseq_100 |>
  rename(External_id = Externel_id) |>
  select(Score, GeneSymbol, External_id) |>
  left_join(cistr_hs_chipseq_meta)

cistrome_hs_chipseq_100 <- cistrome_hs_chipseq_100 |>
  mutate(Factor = str_to_upper(Factor))

cistrome_hs_chipseq_100 |>
  filter(GeneSymbol == 'CCR2') |>
  slice_max(Score, by = c(Factor, Ontology), with_ties = F) |>
  filter(!str_detect(Ontology, 'cell|acute|None')) |>
  mutate(Ontology = str_remove(Ontology, ' \\[.+'),
         cell_factor = str_c(Ontology, '_', Factor),
         cell_factor = fct_reorder(cell_factor, Score, max)) |>
  slice_max(Score, n = 10) |>
  ggplot(aes(Score, cell_factor, color = Ontology, fill = Ontology)) +
  geom_point(size = 2) +
  geom_col(width = .1, position = 'dodge2') +
  scale_color_brewer(palette = 'Paired') +
  scale_fill_brewer(palette = 'Paired') +
  scale_y_discrete(labels = \(x)str_remove(x, '.+_')) +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  labs(title = 'TF with regulatory potential on CCR2 gene',
       subtitle = 'Cistrome ChIP-seq atlas', x = 'ChIP-seq binding score',
       y = 'Transcription factor') +
  theme_jpub

publish_source_plot('ccr2.chipseq.tf')

cistr_hs_chipseq_meta |>
  filter(Ontology == 'None')

cistrome_hs_chipseq_100 |>
  write_csv('00_util_scripts/ref/cistrome_hs_chipseq_beta100.csv')
